﻿/*
Adapted from AAPlus library by PJ Naughter
Permission was granted to convert the library to C#

This source is modified from PJ Naugters AA+ lib
(http://www.naughter.com/aa.html)
and therefore is covered under his copyright license, which
is as follows:

Purpose: Implementation for the algorithms for a Near parabolic orbit
Created: PJN / 21-11-2006
History: PJN / 16-03-2009 1. Fixed a bug in CAANearParabolic::Calculate(double JD, 
                          const NearParabolicObjectElements& elements) in the calculation of the 
                          heliocentric rectangular ecliptical, the heliocentric ecliptical latitude and 
                          the heliocentric ecliptical longitude coordinates. The code incorrectly used the 
                          value "omega" instead of "w" in its calculation of the value "u". Unfortunately 
                          there is no worked examples in Jean Meeus's book for these particular values, 
                          hence resulting in my coding errors. Thanks to Jay Borseth for reporting this bug.

Copyright (c) 2006 - 2009 by PJ Naughter (Web: www.naughter.com, Email: pjna@naughter.com)

All rights reserved.

Copyright / Usage Details:

You are allowed to include the source code in any product (commercial, shareware, freeware or otherwise) 
when your product is released in binary form. You are allowed to modify the source code in any way you want 
except you cannot modify the copyright details at the top of each module. If you want to distribute source 
code with your application, then you are only allowed to distribute versions released by the author. This is 
to maintain a single distribution point for the source code. 

*/

using System;

namespace AAPlus
{
    public class NearParabolicObjectDetails
    {
        public Coordinate3D HeliocentricRectangularEquatorial = new Coordinate3D();
        public Coordinate3D HeliocentricRectangularEcliptical = new Coordinate3D();
        public double HeliocentricEclipticLongitude;
        public double HeliocentricEclipticLatitude;
        public double TrueGeocentricRA;
        public double TrueGeocentricDeclination;
        public double TrueGeocentricDistance;
        public double TrueGeocentricLightTime;
        public double AstrometricGeocenticRA;
        public double AstrometricGeocentricDeclination;
        public double AstrometricGeocentricDistance;
        public double AstrometricGeocentricLightTime;
        public double Elongation;
        public double PhaseAngle;

        public NearParabolicObjectDetails()
        {

        }

        public NearParabolicObjectDetails(
            Coordinate3D HeliocentricRectangularEquatorial,
            Coordinate3D HeliocentricRectangularEcliptical,
            double HeliocentricEclipticLongitude,
            double HeliocentricEclipticLatitude,
            double TrueGeocentricRA,
            double TrueGeocentricDeclination,
            double TrueGeocentricDistance,
            double TrueGeocentricLightTime,
            double AstrometricGeocenticRA,
            double AstrometricGeocentricDeclination,
            double AstrometricGeocentricDistance,
            double AstrometricGeocentricLightTime,
            double Elongation,
            double PhaseAngle
            )
        {

            this.HeliocentricRectangularEquatorial = HeliocentricRectangularEquatorial;
            this.HeliocentricRectangularEcliptical = HeliocentricRectangularEcliptical;
            this.HeliocentricEclipticLongitude = HeliocentricEclipticLongitude;
            this.HeliocentricEclipticLatitude = HeliocentricEclipticLatitude;
            this.TrueGeocentricRA = TrueGeocentricRA;
            this.TrueGeocentricDeclination = TrueGeocentricDeclination;
            this.TrueGeocentricDistance = TrueGeocentricDistance;
            this.TrueGeocentricLightTime = TrueGeocentricLightTime;
            this.AstrometricGeocenticRA = AstrometricGeocenticRA;
            this.AstrometricGeocentricDeclination = AstrometricGeocentricDeclination;
            this.AstrometricGeocentricDistance = AstrometricGeocentricDistance;
            this.AstrometricGeocentricLightTime = AstrometricGeocentricLightTime;
            this.Elongation = Elongation;
            this.PhaseAngle = PhaseAngle;
        }
    }

    public class NearParabolicObjectElements
    {

        public NearParabolicObjectElements()
        {

        }

        public NearParabolicObjectElements(double q, double i, double w, double omega, double JDEquinox, double T, double e)
        {
            this.q = q;
            this.i = i;
            this.w = w;
            this.omega = omega;
            this.JDEquinox = JDEquinox;
            this.T = T;
            this.e = e;

        }

        public double q;
        public double i;
        public double w;
        public double omega;
        public double JDEquinox;
        public double T;
        public double e;
    }

    public static class NearParabolic
    {

        public static double cbrt(double x)
        {
            return Math.Pow(x, 1.0 / 3);
        }

        public static void CalulateTrueAnnomalyAndRadius(double JD, ref NearParabolicObjectElements elements, out double v, out double r)
        {
            double k = 0.01720209895;
            double a = 0.75 * (JD - elements.T) * k * Math.Sqrt((1 + elements.e) / (elements.q * elements.q * elements.q));
            double b = Math.Sqrt(1 + a * a);
            double W = cbrt(b + a) - cbrt(b - a);
            double W2 = W * W;
            double W4 = W2 * W2;
            double f = (1 - elements.e) / (1 + elements.e);
            double a1 = (2.0 / 3) + (0.4) * W2;
            double a2 = (7.0 / 5) + (33.0 / 35) * W2 + (37.0 / 175) * W4;
            double a3 = W2 * ((432.0 / 175) + (956.0 / 1125) * W2 + (84.0 / 1575) * W4);
            double C = W2 / (1 + W2);
            double g = f * C * C;
            double w = W * (1 + f * C * (a1 + a2 * g + a3 * g * g));
            double w2 = w * w;
            v = 2 * Math.Atan(w);
            r = elements.q * (1 + w2) / (1 + w2 * f);
        }

        public static NearParabolicObjectDetails Calculate(double JD, ref NearParabolicObjectElements elements)
        {
            double Epsilon = Nutation.MeanObliquityOfEcliptic(elements.JDEquinox);

            double JD0 = JD;

            //What will be the return value
            NearParabolicObjectDetails details = new NearParabolicObjectDetails();

            Epsilon = CoordinateTransformation.DegreesToRadians(Epsilon);
            double omega = CoordinateTransformation.DegreesToRadians(elements.omega);
            double w = CoordinateTransformation.DegreesToRadians(elements.w);
            double i = CoordinateTransformation.DegreesToRadians(elements.i);

            double sinEpsilon = Math.Sin(Epsilon);
            double cosEpsilon = Math.Cos(Epsilon);
            double sinOmega = Math.Sin(omega);
            double cosOmega = Math.Cos(omega);
            double cosi = Math.Cos(i);
            double sini = Math.Sin(i);

            double F = cosOmega;
            double G = sinOmega * cosEpsilon;
            double H = sinOmega * sinEpsilon;
            double P = -sinOmega * cosi;
            double Q = cosOmega * cosi * cosEpsilon - sini * sinEpsilon;
            double R = cosOmega * cosi * sinEpsilon + sini * cosEpsilon;
            double a = Math.Sqrt(F * F + P * P);
            double b = Math.Sqrt(G * G + Q * Q);
            double c = Math.Sqrt(H * H + R * R);
            double A = Math.Atan2(F, P);
            double B = Math.Atan2(G, Q);
            double C = Math.Atan2(H, R);

            Coordinate3D SunCoord = Sun.EquatorialRectangularCoordinatesAnyEquinox(JD, elements.JDEquinox);

            for (int j = 0; j < 2; j++)
            {
                double v;
                double r;
                CalulateTrueAnnomalyAndRadius(JD0, ref elements, out v, out r);

                double x = r * a * Math.Sin(A + w + v);
                double y = r * b * Math.Sin(B + w + v);
                double z = r * c * Math.Sin(C + w + v);

                if (j == 0)
                {
                    details.HeliocentricRectangularEquatorial.X = x;
                    details.HeliocentricRectangularEquatorial.Y = y;
                    details.HeliocentricRectangularEquatorial.Z = z;

                    //Calculate the heliocentric ecliptic coordinates also
                    double u = w + v;
                    double cosu = Math.Cos(u);
                    double sinu = Math.Sin(u);

                    details.HeliocentricRectangularEcliptical.X = r * (cosOmega * cosu - sinOmega * sinu * cosi);
                    details.HeliocentricRectangularEcliptical.Y = r * (sinOmega * cosu + cosOmega * sinu * cosi);
                    details.HeliocentricRectangularEcliptical.Z = r * sini * sinu;

                    details.HeliocentricEclipticLongitude = Math.Atan2(y, x);
                    details.HeliocentricEclipticLongitude = CoordinateTransformation.MapTo0To24Range(CoordinateTransformation.RadiansToDegrees(details.HeliocentricEclipticLongitude) / 15);
                    details.HeliocentricEclipticLatitude = Math.Asin(z / r);
                    details.HeliocentricEclipticLatitude = CoordinateTransformation.RadiansToDegrees(details.HeliocentricEclipticLatitude);
                }

                double psi = SunCoord.X + x;
                double nu = SunCoord.Y + y;
                double sigma = SunCoord.Z + z;

                double Alpha = Math.Atan2(nu, psi);
                Alpha = CoordinateTransformation.RadiansToDegrees(Alpha);
                double Delta = Math.Atan2(sigma, Math.Sqrt(psi * psi + nu * nu));
                Delta = CoordinateTransformation.RadiansToDegrees(Delta);
                double Distance = Math.Sqrt(psi * psi + nu * nu + sigma * sigma);

                if (j == 0)
                {
                    details.TrueGeocentricRA = CoordinateTransformation.MapTo0To24Range(Alpha / 15);
                    details.TrueGeocentricDeclination = Delta;
                    details.TrueGeocentricDistance = Distance;
                    details.TrueGeocentricLightTime = Elliptical.DistanceToLightTime(Distance);
                }
                else
                {
                    details.AstrometricGeocenticRA = CoordinateTransformation.MapTo0To24Range(Alpha / 15);
                    details.AstrometricGeocentricDeclination = Delta;
                    details.AstrometricGeocentricDistance = Distance;
                    details.AstrometricGeocentricLightTime = Elliptical.DistanceToLightTime(Distance);

                    double RES = Math.Sqrt(SunCoord.X * SunCoord.X + SunCoord.Y * SunCoord.Y + SunCoord.Z * SunCoord.Z);

                    details.Elongation = CoordinateTransformation.RadiansToDegrees(Math.Acos((RES * RES + Distance * Distance - r * r) / (2 * RES * Distance)));
                    details.PhaseAngle = CoordinateTransformation.RadiansToDegrees(Math.Acos((r * r + Distance * Distance - RES * RES) / (2 * r * Distance)));
                }

                if (j == 0) //Prepare for the next loop around
                    JD0 = JD - details.TrueGeocentricLightTime;
            }

            return details;
        }

    }
}
